简捷不简单

导语

北京冬奥会在2025年2月4日开幕，开幕式上雪花的形象令人印象深刻。雪花生长过程能否用复杂系统建模软件来实现？本文用NetLogo软件的元胞自动机模型做了尝试。

研究领域：晶体生长，元胞自动机，多主体建模

纪洪波| 作者

邓一雪| 编辑

1. 问题的提出

从在天安门广场中心布置的“精彩冬奥”主题花坛到冬奥会开幕式演绎一朵雪花的故事[1][2]，北京冬奥会“雪花”元素分外引人注目。

图片来源：https://www.163.com/dy/article/GULGK13205346RC6.html

图片来源：https://www.163.com/dy/article/GVDN6BUH0534A4RJ.html

图片来源：S.Wolfram.A New Kind of Science.Wolfram Media,2002.

推荐书目《A New Kind of Science》中用元胞自动机按照“从一个六边形的黑色元胞开始，如果一个元胞相邻的元胞有黑色的话，这个元胞就变为黑色。”这一简单的规则得到31步斑图。

1月5日仇玮祎老师在自生成结构读书会交流群中推送了《为啥雪花长这样？》一文[3]。再回溯到2018年9月，集智俱乐部发布了《》[4]

在《Netlogo多主体建模》课程中，张江教授主张：对复杂科学最好的入门手段就是自己亲手在计算机中搭建一个复杂系统。Netlogo就是一个非常好的入门工具，它可以让你通过简单的设置和代码编写就能搭建出一个超酷炫的多主体仿真（模拟）程序。

基于以上背景，开始追问“模拟雪花生长的奠基性模型是什么？如何在NetLogo平台上实现仿真？”

2. 模型的梳理

Packard[5]模型的元胞状态为 0 或 1，分别代表气相和固相，网格为三角形网格，每个元胞最近邻都为六个，表现了冰晶胞结构的六角特性。Packard 提出一种常用规则为只有一个邻居元胞为冰相时（记为 n＝1）才能生长，而 n≥2 的情形不能生长；另一个常用规则是用邻居中有奇数个为冰相则生长偶数个为冰相则不生长的规则（即：n＝1、3、5），也称为奇偶规则。虽然Packard模型规则极其简单，却能模拟出基本的板状、柱状和简单六角分枝状的雪花。但由于模型过于简单，其物理意义太过粗略。Gravner & Griffeath（2006）[6]对Packard模型进行了数学分析。

Libbrecht（2005）[7]着重强调，冰的表面层是“准液体”层（quasi-liquid layer），这是解开雪晶生长之谜的主要因素。Gravner & Griffeath（2008）[8]提出的二维模型采用了三种基本的物理量：雪花晶体部分，用a(x，t)表示是否属于晶体，c(x，t) 表示晶体的质量密度；准液体（quasi-liquid），用密度量b(x，t)表示；扩散水汽（vapor），用密度量 d(x，t)表示。在模型中这些量的演化过程可分为 5 类：扩散、凝结、附着、熔化、噪声。二维模型中随着初始密度ρ增加，在适当设置其它参数的前提下，雪晶生长形态表现出的不同特征。

图片来源：

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.121.6149&rep=rep1&type=pdf

图片来源：https://demianrepucci.com/2009/12/31/snowflake-algorithm/

Gravner & Griffeath（2008）[9]基于水汽的扩散、各向异性附着和准液层，引入了一个三维的、便于计算的、介观的雪晶生长模型。研究了若干能对大多数观察到的雪晶体形态忠实地模拟复现的案例。特别是，具有刻面和分支特征的物理标本最引人注目，模型解释了这一现象。模型还复制了许多其他观察到的特征，包括脊、肋、夹层板和空心柱，以及动态不稳定性。观测到的现象一致表明，模型中的因素在物理雪晶体的生长中是至关重要的。上图列出了主要的方程。

模型用C语言编写基本的计算引擎，用MATLAB进行绘图和可视化，用光线跟踪软件POV-Ray进行渲染。下图左边是用Matlab生成的图像，右边是在左图基础上用光线跟踪软件POV-Ray渲染后的效果。

图片来源：

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.136.1150&rep=rep1&type=pdf

3. 程序的设计

Packard模型的元胞状态、三角形网格、六个近邻表现了冰晶胞结构的六角特性，为后续模型所继承。Packard模型奠定了元胞自动机雪花生长模型的基础框架。Wolfram 在《A New Kind of Science》一书中和刘式达，张自国（1989）[10]运用Packard模型进行模拟仿真，为我们编程仿真提供了参考模式。

为构造六个邻居邻域的基本单元，在Netlogo正方形嵌块的基础上用循环嵌套画出有效嵌块组成的网格图。

其基本单元的相对坐标为：

以奇偶规则为例，要实现邻居中有奇数个为冰相则生长，偶数个为冰相则不生长的规则，需区别(t-1)时刻和t时刻状态。为此，在第三象限设置第一象限反对称镜像，按邻域关系在第一象限读取(t-1)时刻周边六个节点嵌块的颜色状态，并计算嵌块颜色为白色的数量作为判据。根据判据在第三象限反对称镜像改变t时刻由海龟坐标确定的当前嵌块的状态。然后用第三象限反对称镜像更新第一象限状态，为下一时刻规则执行做准备。

此外，要注意虚拟空间设置与界面输入必须匹配。

界面

设置

解决了以上三个关键问题后，就可以运用Netlogo四方格相关指令进行编程。程序代码为：

globals [ i j ]

to setup；；初始化

clear-all

create-turtles 1 [ setxy pxcor pycor ]

ask patches

ifelse ( pxcor = 32 and pycor = 32 ) [set pcolor white ] [set pcolor black ]

；；从第一象限copy 到第三象限

ask turtles

let i1 0

while [ i1